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Abstract 

We show how two-dimensional mixed finite element methods that satisfy the conditions of finite 
element exterior calculus can be used for the horizontal discretisation of dynamical cores for nu- 
merical weather prediction on pseudo-uniform grids. This family of mixed finite element methods 
can be thought of in the numerical weather prediction context as a generalisation of the popular 
polygonal C-grid finite difference methods. There are a few major advantages: the mixed finite 
element methods do not require an orthogonal grid, and they allow a degree of flexibility that 
can be exploited to ensure an appropriate ratio between the velocity and pressure degrees of free- 
dom so as to avoid spurious mode branches in the numerical dispersion relation. These methods 
preserve several properties of the C-grid method when applied to linear barotropic wave propaga- 
tion, namely: a) energy conservation, b) mass conservation, c) no spurious pressure modes, and 
d) steady geostrophic modes on the /-plane. We explain how these properties are preserved, and 
describe two examples that can be used on pseudo-uniform grids: the recently-developed modi- 
fled RTO-QO element pair on quadrilaterals and the BOFMl-Pl^jG element pair on triangles. All 
of these mixed flnite element methods have an exact 2:1 ratio of velocity degrees of freedom to 
pressure degrees of freedom. Finally we illustrate the properties with some numerical examples. 

Keywords: Mixed flnite elements, stability, steady geostrophic states, geophysical fluid 
dynamics, numerical weather prediction 
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1. Introduction 



There are a number of groups that have been developing dynamical cores for numerical weather 
prediction (NWP) and climate modelling, based on triangular meshes on the sphere or on the 



dual meshes composed of hexagons together with twelve pentagons (Ringler et al. , 2000 Majewski 



et al. 


2002 


Satoh et al. 


2008 



edge lengths h that satisfy CqH < h < Cih, a.s h 0, where h is the average edge length, for some 
positive constants Cq, Ci. The principal reason for adopting these grids is that they provide a direct 
addressing data structure whilst avoiding the polar singularity of the latitude-longitude grid, which 
introduces a bottleneck to scaling on massively parallel architectures due to the convergence of 
meridians. One approach to developing numerical discretisations on triangular or hexagonal grids 



is to adapt the staggered Arakawa C-grid flnite difference method on quadrilaterals (Arakawa and 



Lamb 


1977) (used in 


Model ( 


Davies et al. 



Lamb, 1977) (used in several currently operational NWP models, such as the UK Met Office Unifled 



functions on the pressure grid that have zero numerical gradient). By deflning discrete curl and 
divergence operators which satisfy div curl= 0, it is possible to construct C-grid discretisations for 



horizontal wave propagation which have stationary geostrophic modes on the /-plane (Thuburn 



et al. , 2009), a necessary condition for accurate representation of geostrophic adjustment processes. 
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These operators can be used to construct energy and enstrophy C-grid discretisations for the 



nonhnear rotating shallow- water equations using the vector invariant form (Ringler et al. , 2010). 



The drawback with using the C-grid finite difference method on triangles or hexagons instead of 
quadrilaterals is that the ratio of velocity and pressure degrees of freedom (DOF) is altered. The 
quadrilateral C-grid has one pressure DOF stored at the centre of each grid cell, and two velocity 
DOF per grid cell (normal velocity is stored at each of the four edges, which are each shared with 
the neighbouring cell on the other side of the face)[^ This is considered the ideal ratio, since the 
velocity then has an equal number of rotational and divergent DOF which are coupled together in 
the correct way so that there are two inertia-gravity modes (the inward and outward propagating 
modes) for each Rossby mode. On the other hand, the triangular C-grid has only 3/2 velocity 
DOF per grid cell, and the hexagonal C-grid has 3 velocity DOF per grid cell. This means that 
the triangular C-grid has four inertia-gravity modes per Rossby mode; the extra spurious inertia- 
gravity branch has a frequency range that decreases with Rossby deformation radius, leading to 
"checkerboard patterns" in the divergence when the deformation radius is small (as it can be in 
the ocean, or when there are many vertical layers). The hexagonal C-grid has an equal number of 
inertia-gravity and Rossby modes; the extra spurious Rossby mode has very low frequencies and 



propagates Eastwards on the /3-plane (Thuburn, 2008). The effects of these spurious Rossby modes 



has not been reported in practice but there are concerns amongst the operational NWP community 
that if spurious modes are supported by the grid, then they might be initialised during the data 
assimilation process or by physics parameterisations (Staniforth, personal communication). It 



may also be the case that the spurious modes lead to spurious spread/lack of spread in ensemble 
forecasts. Careful numerical experiments are required to investigate this concern. 

The finite element method provides the opportunity to alter the number of degrees of freedom 
per triangular element to ameliorate this problem. A number of finite element pairs on triangles 
have been proposed for geophysical fluid dynamics, mostly in the ocean modelling community 



(Walters and CasuUi 



element pair (Brezzi et al. 



1998 



Le Roux et al. 2007). In (Rostand and Le Roux, 2008), the lowest order Brezzi-Douglas-Marini 



Le Roux et al. 1998, 2005 Cotter et al. 2009; Comblen et al. 2010 



1985), known as BDMl, was investigated in the context of the discrete 



shallow-water equations. The velocity space is piecewise linear with continuous normal compo- 
nents, and the pressure space is piecewise constant. The natural data structure for the velocity 
space stores two normal velocity components on each edge, and hence there are 3 velocity DOF 
per triangular element and 1 pressure DOF. There are too many velocity DOF and hence there 
will be too many Rossby modes per inertia-gravity mode, just as for the hexagonal C-grid. 

The key result of this paper is in showing that discretisations of the linear rotating shallow 
water equations on the /-plane constructed using these spaces on arbitrary meshes satisfy a crucial 
property, namely that geostrophic modes are exactly steady. This is achieved by making use of 



the discrete Helmholtz decomposition, within the framework of discrete exterior calculus (Arnold 



et al. , 2006). As described in (Arnold, 2002), existence of such a decomposition requires that the 



following diagram commutes: 



(1) 



E 



S 



V 



^Here, and in the rest of the paper, we consider compact domains without boundary such as the sphere and 
rectangles with double periodic boundary conditions. 



2 



where 11^, 11'^ and 11^ are suitably chosen projection operators. The same Helmholtz decompo- 
sition can then be used to study the discrete dispersion relations for the numerical discretisation. 
Within this framework, we then conclude that an optimal choice is to have dim(5') = 2 dim(l^) 
which, at least in the periodic plane, satisfies necessary conditions for absence of both spurious 
inertia-gravity and spurious Rossby waves. 

The rest of this paper is organised as follows. The general framework of mixed finite element 
methods applied to the linear rotating shallow-water equations is described in Section [2} and the 
four properties of energy conservation, local mass conservation, absence of spurious pressure modes 
and steady geostrophic modes are discussed. In Section |3| two examples are then introduced that 
fit into this framework, and numerical results are presented in section |4j Finally, we give a summary 
and outlook in Section IH 



2. Mixed finite elements for geophysical fluid dynamics 

In this section we describe how mixed finite elements can be used to build flexible discretisa- 
tions on pseudo-uniform grids. We concentrate on the rotating shallow-water equations which are 
regarded in the numerical weather prediction community as being a simplified model that contains 
many of the issues arising in the horizontal discretisation for dynamical cores. Since in this paper 
we are concerned with wave propagation properties, we restrict attention to the linearised equations 
on the /-plane, /3-plane or the sphere. First, we introduce the mixed finite element formulation 
applied to the linear rotating shallow-water equations, then we discuss various properties of the 
formulation that are a requirement for numerical weather prediction applications, namely global 
energy and local mass conservation, absence of spurious pressure modes and steady geostrophic 
states. These properties all rely on exact sequence properties, i.e. div-curl relations, as described 



in (Arnold et al. 2006). 



2.1. Spatial discretisation of the linear rotating shallow-water equations 

In this paper we consider the discretisation of the linear rotating shallow-water equations on a 
two dimensional surface fl that is embedded in three dimensions (which we restrict to be compact 
with no boundaries, e.g. the sphere or double periodic x — y plane): 

Ut + fu^ + c^Vrj = 0, r]t + V ■u = 0, u ■ n = on dfl, (2) 

where u = {u, v) is the horizontal velocity, u-^ = k x u, f is the Coriolis parameter, = gH, g is 
the gravitational acceleration, H is the mean layer thickness, h = H{1 + 77) is the layer thickness, 
k is the normal to the surface Q, and V and V- are appropriate invariant gradient and divergence 
operators defined on the surface. We form the finite element approximation by multiplying by 
time-independent test functions w and 0, integrating over the domain, integrating the pressure 
gradient term c^Vrj by parts in the momentum equation, and finally restricting the velocity trial 
and test functions u and it> to a finite element subspace S C H{dw) (where if(div) is the space 
of square integrable velocity fields whose divergence is also square integrable), and the elevation 
trial and test functions rj and a to the finite element subspace V G L'^ (where is the space of 
square integrable functions): 



d 

di 



[ w'' ■ dV + [ fw^ ■ (u^)^ dV - [ V ■ w'^T]^ dV = 0, \/w''eS, (3) 
Jn Jn Jn 

[ aVdV^+ / a^'V-u^dV = 0, Va^ G V. (4) 
Jn Jn 



d_ 
di 
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After discretisation in time, these equations are solved in practise by introducing basis expansions 
for w^, u^, 1]'^, and a'^ and solving the resulting matrix-vector systems for the basis coefficients. 
In this framework we restrict the choice of finite element spaces 5* and V so that 



u 



h 



eS ^ v-u^ eV. 



The divergence should map from S onto V, so that for all functions (p^ & V there exists a velocity 
field & S with V ■ ■u'^ = (f)^. Such spaces are known as "div- conforming" . Furthermore we 
require that there exists a "streamfunction" space E C such that 

where the A; x V operator (the curl, which we shall write as V"*") maps onto the kernel of V- in 
S. A consequence of these properties is that functions in E are continuous, vector fields in S only 
have continuous normal components and functions in V are discontinuous. 

2.2. Energy conservation 

Global energy conservation for the linearised equations is a requirement of numerical weather 
prediction models for various reasons, in particular because it helps to prevent numerical sources of 
unbalanced fast waves. It is also a precursor to a energy-conserving discretisation of the nonlinear 
equations using the vector-invariant formulation. For the mixed finite element method, global 
energy conservation is an immediate consequence of the Galerkin finite element formulation. The 
conserved energy of equations ^ is 

H = - I |wp + cVdV. 
2 Jn 

Substituting the solutions and rj^ to equations and taking the time derivative gives 

d 



-H=[ ■ ii^ + c^r]''f]^ dV. 
' Jn 



dt 

Choosing = and = rj^ in equations (|3||4|) then gives 
d 



dt 



H = I ■ v!" + C^7]^7]^ dV 

Jn 

= J -fu''- {u^)^ + c^V-u''ri''-c^r]^V-u\ dV = 0. 



=0 =0 

2.3. Local mass conservation 

Local mass conservation is a requirement for numerical weather prediction models since it 
prevents spurious sources and sinks of mass. For the nonlinear density equation, this can be 
achieved using a finite volume or discontinuous Galerkin method. For mixed finite element methods 
of the type used in this paper applied to the linear equations, consistency and discontinuity of 
functions in V requires that element indicator functions {i.e. functions that are equal to 1 in 
one element and in the others) are contained in V. Selecting the element indicator function for 
element e as the test function in equation Q gives 

[ r]''dV+ [ u^-ndS = 0, 

dtJe Jde 

where de is the boundary of element e. Since has continuous normal components on element 
boundaries, this means that the flux of rj'^ is continuous and hence rj'^ is locally conserved. 
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2.4- Absence of spurious pressure modes and stability of discrete Poisson equation 

The principle reason for using the staggered C-grid for numerical weather prediction is that 
the collocated A-grid, in which pressure and both components of velocity are stored at the same 
grid locations, suffers from a checkerboard pressure mode which has vanishing numerical gradient 
when the centred difference approximation is used, despite being oscillatory in space. This pressure 
mode rapidly pollutes the numerical solution in the presence of nonlinearity, boundary conditions 
and forcing, and can be easily excited by physics subgrid parameterisations or initialisation via 
data assimilation from noisy data. 

In the context of mixed finite element methods applied to the equation set (|2]), spurious pressure 
modes relate to the discretised gradient Dcj)^ G S' of a function (p^ G V defined by 



w'' ■ D<p^ dV = - I V ■ w^<p^ d V, \/w^ G S. 
'n 



On uniform grids, spurious pressure modes are functions (p^ from the pressure space V which 
have zero discretised gradient D(j)^ even though V(f)^ is non-zero. On unstructured grids or grids 
with varying edge lengths, spurious pressure modes are functions which have discretised gradient 
becoming arbitrarily small as the maximum edge length ho tends to zero, despite their actual 
gradient staying bounded away from zero. Such functions would prevent the numerical solution 
of equations ^ converging at the optimal rate predicted by approximation theory. We make the 
following definition of a spurious pressure mode. 

Definition 1 (Spurious pressure modes). A mixed finite element space {S,V) is said to be free of 
spurious pressure modes if there exists 72 > independent of ho such that for all (j)^ G V , there 
exists nonzero E S satisfying 

1 0"V-i;'^dF>72||0'^|U,||«"||^^(div). (5) 
Jo, 

Condition ^ is one of two sufficient conditions for numerical stability of the mixed finite 
element discretisation of the Poisson equation — = / given by 



-v^dV = - / V ■ w^^^ d V, \fw^ G S, 
n Jq 

Jn Jq 

This discretisation is stable [i.e. small changes in the right-hand side lead to small changes in the 
solution field in the limit as ho tends to zero) if Condition ^ holds, together with the condition 
that there exists 71 > independent of ho such that 

/ v'^-v'^dx>^,\\vYHi,^), (6) 
Jn 



for all e S such that /V ■ v'^(f)'^dV = for all (f)'^ G V. As reviewed in Arnold (2002), 
Condition (|5| is satisfied if it is possible to define a bounded projection 11'^ : H(div) — )■ S such 
that the following diagram commutes: 



H{div,n) L2{n) 



S V 



(7) 
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where 11^ is the usual L2 projection operator. This means that taking any square integrable 
velocity field u with square integrable divergence, evaluating the divergence and projecting into 
V produces the same result as projecting u into S using 11'^ and evaluating the divergence. The 
projection 11'^ is constructed by applying an projection of normal components on element edges, 
ensuring that u is L^-orthogonal to gradients of functions from V in each element, and ensuring the 
remaining degrees of freedom in u are L^-orthogonal to divergence-free functions in each element. 
We shall explain how this is done for the two examples described in Section [3j To check that the 
diagram ([T]) commutes, it is sufficient to show that 

f a^{V -u-V ■Ti^u)dV = Wa^ eV,u e H{div,K), 

for each element K, since this defines the L2 projection 11^ into the discontinuous space V. This 
is easily checked using integration by parts: 

[ a''V-udV = - [ Va''-udV+ I a^u-ndS, 
Jk Jk JdK 

= - I Va^ ■ li^u dV+ I a^U^u -ndS = [ a^V ■ li^u d V, 
Jk JdK Jk 

as required. 

As also reviewed in Arnold (2002), Condition ^ is satisfied if vector fields v E S with diver- 



gence orthogonal to V are in fact divergence-free. This is satisfied by the types of mixed finite 
element methods considered in this paper since the divergence maps from S into V, and so the 
projection of V ■ into V is simply the inclusion. Hence, if the divergence is orthogonal to V, 
the divergence must be zero, and so ^ is satisfied. 

2.5. Discrete Helmholtz decomposition 

Proof of the condition that geostrophic modes are steady requires the construction of a discrete 
Helmholtz decomposition. Since Condition ^ holds, the discrete gradient operator D : V ^ S, 
has no non-trivial kernel. For any e E, the curl V"*" of '0'* satisfies 



f V^^Ij^ ■ Dcj)^ dV = - I V ■V^i)^(t)^dV = 
Jn Jq^ ^ ' 



for any (p^ E V, and hence the curl from E to S and the discrete divergence from V to S map onto 
orthogonal subspaces of 5*. This means that there is a one-to-one mapping between elements of S 
and E X V, defining a discrete Helmholtz decomposition 



u 



h - ^^^h ^ jj^h ^^h^ ^ ^ ^E,(P^ E V, E H, (8) 



where H G S is the space of discrete harmonic velocity fields 

H'' = !^u^ E S -.V ■u'' = 0, j - V^tlj^ dV = Q.-iilj^ E E^ . 
The dimension of is the same as the dimension of the space H of harmonic velocity fields 
= |u e i/(div) : V = 0, j -V^iljdV = Q,^il) E H^^ , 
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i.e., velocity fields with vanishing divergence and (weak) curl (In the periodic plane, these harmonic 
velocity fields are the constant velocity fields, but there are no harmonic velocity fields on the 



sphere); however ^ H in the general case (Arnold et al. , 2006). The kernel of V"*" in E is 



the subspace of constant functions, and stability results (as described in Section 2.4) imply that 
the kernel of D in ^ is the subspace of constant functions, and hence we can use Equation ^ to 
obtain a DOF count for 5*. 

dim(5) = (dim(E) - 1) + {dimiV) - 1) + dim(iJ), 

and hence 

dim(^) = dim(5) - dim(l^) + 2 - dim(i7). 
For our DOF requirement dim(5') = 2dim(\/), we obtain 

dim(E) = dim(r) + 2 - dim(/f), 

which becomes dim(£') = dim(y) for the periodic plane and dim(i?) = dim(\^) + 2 for the sphere. 
If dim(S') > 2dim(V^), then dim(i?) > dim(V^) + (2 — dim(if)) and vice versa. This will become 



important when we examine wave propagation in Section 2.8 



2.6. Vorticity and divergence 

The discrete vorticity associated with the velocity u'^ & S is defined as G £" such that 



^^^'^ dv = - v^r ■ w d V, y-f" G E. (9) 

n Jn 

It is possible to obtain u E S from the discrete vorticity C, E E and the divergence 6^ = V -u^ E V 
by solving two elliptic problems for the streamfunction ip^ and velocity potential (p^. To obtain 
the streamfunction G E, we use the Helmholtz decomposition and rewrite equation (|9| as 

n 





/ 




Jn 



I ^'^dy = o, 

Jn 



which is the usual finite element discretisation of the Poisson equation for tp^. To obtain the vector 
potential 0^ requires the solution of the coupled system 

[a^V-D(j)^dV = [ a^S^dV, Wa^ e \a:aeV, [adV = o\, 
Jn Jn \^ Jn J 

[ w'' ■ D(l)'' dV = - f V ■w^cjl'dV, yw'^ES, [(l)^dV = 0. 
Jn Jn Jn 

This is the mixed finite element approximation to the Poisson equation already discussed in Section 
2.4 if the Conditions (|5| and ^ are satisfied, the coupled system is well-posed. 

2. 7. Steady geostrophic modes 

On the /-plane (planar domain with constant /), geostrophic balanced states satisfying fw^ + 
c^Vt] = are steady since V ■ u = 0. The remaining solutions of the linear rotating shallow- 
water equations are fast inertia-gravity waves. In the quasi-geostrophic limit (slow, large scale 
motion), when nonlinear terms and spatially varying / are introduced, these steady states become 
slowly-evolving balanced states that characterise large-scale weather systems. It is crucial that 
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a discretisation gives rise to steady geostrophic states on the /-plane, otherwise when nonhnear 
terms and spherical geometry are introduced, balanced states will emit noisy inertia-gravity waves 
that will pollute the numerical solution over timescales that are much shorter than that required 
for a weather forecast. To show that mixed finite element methods have steady geostrophic modes. 



we follow the approach of Thuburn et al. (2009 ), namely we aim to show that vanishing divergence 
implies steady vorticity, then checking that vanishing divergence and steady vorticity implies steady 
velocity. 

To obtain a geostrophic balanced state corresponding to a given streamfunction ip'^, we initialise 
u'^ and 7]^ as follows: 

1. Set u'' = V^ip^. 

2. Set rf^ from the geostrophic balance relation 

/ aVdV^ = / / aVdr, Va^ G V. (10) 
Jo, Jn 

Substitution in equation ^ then gives 



d 

di 



n Jn Jn 

= f [ V-ioV^dV-c^ [ V-i^VdV, 
Jn Jn 

= 0, 



having noted that V ■ G V and so we may choose a'^ = V ■ in equation (10). To show that 



f]'^ = 0, first note that u'^ = V-^il)^ and hence V • w'^ = 0. Equation ^ thus becomes 

a^i]'' d = 0, Ma^ e V, 



n 



and hence r/^ = 0. This means that the geostrophic balanced state is steady. 

2.8. Numerical dispersion relations 

In this section we consider the numerical wave propagation properties of this family of finite 
element discretisations, on the /-plane and on the /3-plane in the quasi-geostrophic limit. 

Dispersion relations are computed by assuming time-harmonic solutions proportional to e"*"^* (a 
valid assumption if the equations are invariant under time translations) and studying the resulting 
eigenvalue problem. For the continuous equations on the periodic plane, the equations are also 
invariant under spatial translations and so it may be assumed that the eigensolutions take the form 
^i(k-x-ujt) ^j-^gj-g is restricted so that the periodic boundary conditions are satisfied. Substitution 
in the equations of motion leads to an algebraic system relating k to u: the dispersion relation. 
For the linear shallow-water equations this system is most easily obtained by using the Helmholtz 
decomposition for u. Numerical dispersion relations for continuous-time spatial discretisations are 
also computed by assuming time-harmonic solutions, leading to a discrete eigenvalue problem. If 
a structured mesh is used on the periodic plane with a set of discrete translation symmetries then 
eigensolutions satisfy the property that translating from one cell to another by Ax results in the 
discrete eigensolution changing by a factor of e*^*'^^^ where k is again chosen so that the periodic 
boundary conditions are satisfied. This can again lead to a numerical relationship between k and 
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u, obtained for both the /-plane, and the /3-plane in the quasi-geostrophic hmit, for the hexagonal 



C-grid in Thuburn (2008), and for the PIdg — -P2 finite element pair in Cotter and Ham (2011). 



Here, we discuss the properties of the discrete eigenvalue problem arising from the finite element 
spaces from the framework of this paper. The discussion makes use of the discrete Helmholtz 
decomposition. In the /-plane case, substitution of the discrete Helmholtz decomposition into 
equations (pH and assuming time-harmonic solutions yields 



-ico [ V-f'' ■V^j'^dV + [ fV^^-D(t)^dV = 0, (11) 
Jn Jn 

-icu [ Da^-D(P^dV+ [ fDa^-h^^+{D(l)^y)dV-c^ [ V ■ Da\^ dV = 0, (12) 
Jn Jn ^ ' Jn 

-iu [ aVdV"+ / a^V-D^^dV = 0, (13) 
Jn Jn 



for all test functions a G V, 7 G E. Next we define projections P : V ^ E and P : E ^ V 
by 

I V^^ -V {P^<p^)dV = [ V-f^ ■ D4>^ dV, W4>^eV, -f'^eE, 
Jn Jn 

[ Da^ ■ D{P^ij^)dV = [Da^-V^^dV, e E, e V. 

Jn Jn 

These projections are uniquely defined since P^ uses the standard continuous finite element dis- 
cretisation of the Laplace operator which is solvable by the Lax-Milgram theorem when E is re- 
stricted to mean zero functions, and P^ uses the mixed finite element discretisation of the Laplace 
operator using the spaces S and V which is solvable by the stability conditions ^ and ^ when 
V is also restricted to mean zero functions. 

Using these projections, and the fact that the divergence operator maps from S to V, equations 



(11 13) become 



-iuj^lj^ + fp^<i)^' = 0, 

-iu I Da^-D(p^dV + f [ Da^ ■ DP^ij^dV 
Jn Jn 

+ [ fDa^ ■ {Dcj)^)^ dV -c^ f V-DaVdV" = 0, 
Jn Jn 

-iujT]^ + V ■ D(p^ = 0, 
and elimination of ip^ and use of the definition of D gives 

= uj({uj^ + f) [ aVdV^+ / aVdV^-c^ f V ■ Da^r]^ dV 
\ Jn Jn Jn 

+if I Da^ ■D{P^P^(t)^ -(l)^)dV -Lo I fDa^ ■ {Dcj)^)^ dV, 
Jn Jn 



(14) 

(15) 
(16) 



(17) 



where (j)^ is obtained from equation (16). The first row of equation (17) is the discretisation of 



the continuous eigenvalue problem for the rotating shallow-water equations using the mixed finite 
element spaces V and S. In this case the eigenvalues of this discrete eigenvalue problem converge 



to the eigenvalues of the continuous problem at the optimal rate as described in Boffi et al. (1997). 
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However, there are two extra terms in the bottom row of equation (17). The second term converges 
to zero for smooth 0'^, and use might be made of spectral perturbation theory to examine what 
effect this has on the discrete eigenvalue problem; we have not yet developed a technique to do 
this. However, the impact of the first term in the second row is more immediately clear, since 
it involves projecting (p^ from ^ to ii^ and back to V again. If V has larger dimension than E, 
which is the case for the lowest order Raviart- Thomas element on triangles, for example, then this 
double projection will have a kernel, and (p^p^ — 1)0'^ will not be small. This leads to spurious 
branches of inertia-gravity waves, i.e. branches of solutions of the discrete eigenvalue problem that 



do not converge to solutions of the continuous eigenvalue problem as /i — )■ 0. See Danilov (2010) 
for numerical examples illustrating this spurious modes, in particular Figures 2,3 and 4. Hence, 
dim(y) < dim(£') is a necessary condition for the absence of spurious divergent inertia-gravity 
modes. 

A similar approach can be taken to studying the /3-plane solutions in the quasi-geostrophic 
limit. Substitution of the discrete Helmholtz decomposition into equations (|3]|4]) and assuming 
time-harmonic solutions yields 

-iuj I V-i^ - Vij^dV + I {h + l3y)V-i^ ■ D(t)^dV = 0, (18) 
Jn Jq 

-icu [ Da'' ■ {D(l)^ + V^i;^)dV+ 
Jn 

j if + /3y)Da'' ■ (Vip'' + {D(l)'')^^dV - j V ■ Da''r]''dV = 0, (19) 

-ico [ aVdV"+ / a'^V-D^^'dV = 0. (20) 
Jn Jn 

In the usual quasi-geostrophic limit, the leading order solution is 

0j = 0, [ foDa'' ■ VV-J dV + c^ [ V ■ Z^a'^r/J dV = 0, 
Jn Jn 

where (pg, ipg and r]g are the leading order terms in the low Rossby number expansion of 0'*, ip'^ 
and respectively. This is the same as the geostrophic steady state formula for the /-plane, and 
we have 

The next order in the expansion of the equations (we do not make use of the next order in the 
equation) is 

-zu [ V7"-VV'Jd\/+ / /oV7^-/^0^,d\/+ / /3y V7^ ■ V^^-J dl^ = 0, (21) 
Jn Jn Jn 



-ico [ a\Jdr+ / a'^V ■ Dcp'^^gdV = 0. (22) 
Jn Jn 



Again, the embedding property implies that icorig = V ■ Dcp'^g. Since is continuous and D4>'^g has 
continuous normal components, we may integrate by parts in the second two terms in equation 



(21 ), to obtain 

= -ico [ V-f'' -Vtp^dV -ico [ ^j'^tp'^dV- [ l3-f''^i^^dV 
Jn Jn c Jn ox 



Jn c 



10 



The first line is the continuous finite element approximation to the Rossby wave eigenvalue problem 
using the finite element space E, which has convergent eigenvalues. The second line is a pertur- 
bation involving (l — p^p^^ which will not always be small if P^P^ has a non-trivial kernel. 
This will be the case if dim(l^) < dim{E), as occurs in the lowest order Brezzi-Douglas-Marini 



(BDMl) element on triangles (Brezzi et al. , 1985) which has PI as the streamfunction space, and 
hence 2dim(\/) = dim{E) + 2-dim{H). If P^P^ has a non-trivial kernel, this will lead to spurious 
Rossby wave branches of the numerical dispersion relation. We conclude that dim(V^) = dim{E) is 
a necessary condition for avoiding both spurious divergent modes and spurious irrotational modes. 
Note that this is not a sufficient condition since it is still possible for p^p^ or p^p^ to have 
non-trivial kernel even in this case. This condition motivates the selection of examples of mixed 
finite element spaces given in the next section. 

3. Examples 

In this section we provide two examples of mixed finite element spaces that are suitable for 
constructing pseudo-uniform grids on the sphere, and that have the additional property that there 
are exactly twice as many velocity degrees of freedom as pressure degrees of freedom, which prevents 
the presence of spurious mode branches. The first example is the modified Raviart-Thomas element 
on quadrilaterals, and the second example is the Brezzi-Douglas-Fortin-Marini element on triangles. 

3.1. Modified Raviart-Thomas element on quadrilaterals 

There have been several efforts at developing numerical weather prediction models based on a 



cubed sphere grid (see Putman and Lin (2007), for example) in which a grid on the surface of a 



cube is projected to a sphere. The drawback in using such is grid is that to obtain a C-grid finite 



difference method with stationary geostrophic states, the scheme of Thuburn et al. (2009) must 
be used, which requires the grid to be orthogonal in the sense that lines joining adjacent pressure 
nodes must cross cell boundaries at right-angles. On the cubed sphere, this condition does not 
produce a pseudo-uniform grid since elements become clustered near the poles as the resolution 
is increased. Mixed finite elements provide extra freedom to design numerical schemes since the 
orthogonality condition is not a requirement; it is replaced by the conditions on finite element 
spaces specified in Section [2j 

The lowest-order Raviart-Thomas finite element space is the mixed finite element analogue 
of the C-grid since the pressure space is piecewise constant functions, and the velocity fields are 
constrained to be have constant, continuous normal components on element edge. This means that 
one normal component of velocity must be stored on each element edge, just like the C-grid. The 
velocity fields are constructed on a square 1x1 reference element K with coordinates (^1,^2), on 
which the ^i-component of velocity ii is obtained by linear interpolation between constant values 
on the ^1 = and ^1 = 1 edges, and the ^2-component is obtained by hnear interpolation between 
constant values on the ^2 = and ^2 = 1 edges. In these coordinates, the divergence is constant. In 
any physical element K in the mesh, we define a coordinate mapping gr : ^ 1— > x, and the velocity 
in K is obtained via the Piola transformation 

^(^) = T^lf •^(^)' 

det (I) 9^ 



which preserves fiux integrals 



f u-ndS{$) = I u-ndS{x 

Jj Jgh) 
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guaranteeing continuity of normal fluxes. The divergence satisfies 

where V is the divergence in the local coordinates If the coordinate transformation is affine 
(elements are parallelograms), the determinant of the Jacobian is constant, and so the divergence 
of the velocity is constant in each element. However, for general quadrilateral elements (required 
for the cubed sphere), the coordinate transformation is bilinear, with linear determinant of the 



Jacobian. The solution, proposed by Boffi and Gastaldi (2009), is to modify the basis functions by 



adding a divergent correction with vanishing normal components on the boundary that makes the 
divergence constant. The corresponding streamfunction space E is the usual continuous bilinear 
space on quadrilaterals, often denoted Ql, and it can easily be shown that the V"*" operator maps 
from E into S in this case. In fact, the Bofii-Gastaldi correction adds a purely divergent component 
to the velocity field and so the V"*" embedding property is not affected. 

The RTO-QO finite element space has one pressure degree of freedom per quadrilateral element, 
and one velocity degree of freedom per edge. Since (for periodic boundary conditions or the 
sphere) each edge is shared by two elements, this means that there are exactly twice as many 
velocity degrees of freedom as pressure degrees of freedom. This modified Raviart-Thomas finite 
element space satisfies all the conditions that we require in this paper and hence has potential for 
use on pseudo-uniform grids for numerical weather prediction. 

3.2. Brezzi-Douglas-Fortin-Marini element on triangles 

There is an analogous Raviart-Thomas finite element space on triangles which satisfies the 
required embedding properties. However, these spaces satisfy 2 dim(l^) > dim(S') in general. For 
example, the lowest order finite element space RTO-PO has one pressure degree of freedom per 
element, and one velocity degree of freedom per edge, meaning that 3dim(V^) = 2dim(S'). The 
BDMl element on triangles has one pressure degree of freedom per element and two velocity 
degrees of freedom per edge, meaning that 3dim(\^) = dim(5'), so 2dim(\^) < dim(5'). However, 
the little-used lowest order Brezzi-Douglas-Fortin-Marini (BDFMl) element together with PIdg 
on triangles satisfies 2dim(l^) = dim(S'). The BDFM family of elements for quadrilaterals was 



introduced in Brezzi et al. (1987), and an analogous family for triangles was described in Brezzi 
and Fortin (1991). On triangles it is infrequently used since the BDM and RT families have less 
degrees of freedom for the same order of convergence (after suitable post-processing). However, 
these extra degrees of freedom are useful to us here since they mean that dim(l^) = dim(£'). 

Here we describe the BDFMl element on triangles as an augmentation of the BDMl element 
on triangles, which we recall first. Given a triangle K, we define Pk{K) to be the space of /c-th 
order polynomials on K. We define the following spaces on K: 

velocity space S{K) = {Pi{K)y 
pressure space ^(-^) = Po{K). 

For a triangulation T of the domain Q, we define the BDMl velocity space 

S = {ve H{diY, n):v\Ke S{K), K e T}, 

where H{div, Q) is the space of vector fields with square integrable divergence, which requires that 
V has continuous normal component across triangle edges. The pressure space is 

V = {v:v\KeViK)}, 
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with no continuity requirements across edges. 

A convenient set of local nodal basis functions for S is defined by choosing two node points on 
each triangle edge, each node located at one of the vertices belonging to that edge: a total of six 
node points. For example, in the triangle shown in Figure [T| on edge el there are two node points, 
one at vertex v3 and one at vertex v2. The basis function associated with edge el and vertex v3 is 

01,3 = ^2-^3, 

where ^2 is the unit tangent vector to edge e2 and where {Xi}^^^ are the barycentric coordinates 
associated with vertices el, e2 and e3 respectively. It can easily be checked that 4>i 3 has normal 
component equal to 1 at the node point located at vertex v3 on edge el, and normal component 
equal to zero at each of the other node points. The other six basis functions are constructed in a 
similar manner. 

To increase the number of degrees of freedom in each triangle K in the triangulation T, we 
define the local BDFMl space S{K) by 

S{K) = {v e P2{Kf : V • n = on dK}. 

Since all of the vectors in S{K) vanish on the boundary of they do not alter the values of the 
normal components at the boundary, and so there are no additional continuity constraints. The 
dimension of P2{K)'^ is 12, and there are 9 independent degrees of freedom which do not vanish 
on the boundary, which means that dim(S'(i^')) = 3. 

A convenient set of local nodal basis functions for S is defined by locating nodes that store the 
tangential component of velocity at the centre of each edge. The tangential component of velocity 
is permitted to be discontinuous and so a different value of the tangential component will be stored 
on each side of the edge. The basis function associated with the node at the centre of edge el is 

4>i = 4tiA2A3. 

It can easily be checked that 4>i has vanishing normal component on all edges, tangential component 
equal to 1 at the centre of edge el and vanishing tangential component on the other two edges. 
The other two basis functions are constructed in a similar manner. 
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The augmented velocity space 5* on the triangulation T is defined as 



V 



i), v'\k e SiK), vIkE SiK), K e T}. 



S = {ve H{diY,n) : V 
The pressure space V is defined as 

V = {r]E Pi{K)} 

with no continuity requirements. For this mixed element pair the velocity space S has 6 DOF per 
element, and the pressure space V has 3 DOF per element, hence there are twice as many velocity 
DOF as pressure DOF, just as for the C-grid finite difference method on quadrilaterals. 

For our augmented velocity space, it is easy to define the projection operator 11'^. The projection 
is computed element by element and guarantees the continuity oi u - n across element edges. The 
projection on an element K is defined from the following conditions: 

/ -^^{Tl^u-u) -ndS = y-f'' E P\e{i)),y edgese{i) edK,i = 1,2,3, (23) 

Je{i) 

I V-i^ ■{Tl^u-u)dV = y-f''EP\K), (24) 
Jk 

V^B ■ (n^u-u)dV = 0, (25) 

K 



where B is the cubic "bubble" function (as used in the MINI element (Arnold et al. , 1984[ )). In a 
triangle K, the cubic bubble function B^ is the unique cubic polynomial which takes the value 1 
at the barycentre and on all three edges. The streamfunction space E is 

E = {tPE H\Q) : = ^'\k + aBK, ^'\k EP2{K),aE R}. 



Equation (23) comprises the BDMl projection operator, fixing six degrees of freedom. The com- 
ponents of the extra degrees of freedom S{K) are not affected since they all satisfy u ■ n = on 
dK. The vector field V"*"!? lies inside S{K) since it is quadratic (being the skew gradient of a cubic 
function, B) and has vanishing normal component on dK (since B is zero on dK). If we construct 
an orthogonal (relative to the L2 inner product) decomposition of S{K) into V"*"!? © S{K) then 



we see that equation (25) only involves the V -B component, and equation (24) only involves the 



remaining two S{K) components, as 



f V-f^ -V^BdV = - [ ■V-f''BdV+ [ Vj^-nBdS, 

Jk Jk^ ^ JdK 

=0 =0 

because B vanishes on dK. The space {v = V'')^,'-)^ E P1{K)} is spanned by constant vector 
fields, and hence equation (24) fixes the two degrees of freedom in S{K). Bounds on 11'^ can be 



obtained by following the steps of Brezzi et al. (1985), since it simply involves L2 projection onto 
various moments. 

We define the streamfunction space E as the usual Lagrange continuous quadratic space aug- 
mented by cubic bubble functions. For any function ip E E, the curl V"*" maps into S: V"*"^ E S. 
Furthermore, we may define a projection operator II'^ : H^{Q) — )■ H{div) by 



n^^d^ 

U^tpdV 



'?/'(t>j) V vertices fj, i = 1,2,3, 
ipdS, Vedgescji = 1,2,3, 

V'dV, 



K 



K 
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for each element K. To show that the projections commute with V"*", i.e. U^V'^tjj = V'^H^tjj, we 



check each of the conditions (23 25). Condition (23) becomes 



Je{i) Jeii) 



ijV-f^dx + [-f''ij] 



e(i) 



e(i) ' 



e{i) 



f 7^V^n^^-nd^, y-f^ e P\eii)), i = 1,2,3, 

Je{i) 



(26) 



where v^^-^ are the two vertices at either end of edge e{i), and having noted that is constant 



for 7^ G P^{e{i)). Condition (24) becomes 



V7'' ■ n^v^v^ d V 



K 



Jk 

- [ -f'^V -V^tpdV + [ -f^V^^-ndS, 
Jk ^ ^ ^ JdK 



V7'^- v^n^Vdi^, y-f''eP\K) 



K 



where we have used equation (26). Finally, condition (25) becomes 



K 



K 



K JdK " 



K 



V^B ■V^Ii'^i)dV, 



K 



since V^-B is constant in K and B is zero on dK. 
Counting global degrees of freedom, 

dim(E) = A^edge + N^ert + ATf^ee = 2N^A^^ + C, dim(5) = 3Aredge, dim(V) = 2N, 



edge 



3A^i 



face; 



where C is the Euler characteristic of the domain Vt which is equal to for the doubly-periodic 
domain and equal to 2 on the sphere. On the sphere there are two extra constraints: namely that 
the divergence and the vorticity both integrate to zero, and so in both cases dim(i?) + dim(l^) = 
dim(S'). Finally, we note that each triangle has three edges which are each shared with one other 
triangle, and hence 2Nedge = SiVface- 
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4. Numerical results 



In this section we illustrate the properties of the BDFMl finite element space applied to the 
linear rotating shallow-water equations. The equations were integrated numerically using the 
implicit midpoint rule, and the resulting discrete system was solved by using hybridisation which 



is a standard technique for solving elliptic problems (see Brezzi and Fortin (1991) for a detailed 
description) in which the continuity constraints on the velocity space are dropped, and are instead 
enforced in the equation by Lagrange multipliers. It becomes possible to eliminate both the 
velocity and free surface variables from the matrix equation, leaving a symmetric positive definite 
system to solve for the Lagrange multipliers. The velocity and free surface variables can then be 
reconstructed element-by-element. One of the benefits of this approach is that it can be applied 
when the Coriolis term is present, resulting in a fully implicit treatment of this term. In our 
numerical tests this system was solved using a direct solver. In the case of BDFMI-PI/jg, there 
are three Lagrange multipliers per element. 

In the test cases with variable Coriolis parameter /, a continuous piecewise quadratic repre- 
sentation of / was used. 

4-1- Steady states for the f -plane 

We verified that the geostrophic states are exactly steady on the /-plane for the BDFMl finite 
element space by randomly generating streamfunction fields ip from the streamfunction space S 



on the same mesh as used for the PIdg — -P2 finite element pair steady state tests in Cotter et al. 



(2009), with streamfunction equal to zero on the boundary. This mesh is a planar unstructured 
mesh in the x — y plane in a 1 x 1 square region. The velocity was initialised by setting u = kx Vip 
where k is the unit normal to the domain i.e. k = (0,0, 1), and rj was obtained by solving the 
discrete elliptic system 

w''-v^dV+ [ c^V-icVdV^ = (27) 
n Jn 



[ a''V-vdV = [ Da'' ■ f (u^)^ dV, 
Jq Jq 



{2i 



with = / = 1. We then integrated the equations forward for arbitrary lengths of time and 
observed that the layer thickness h and velocity u remained constant up to machine precision. We 
also conducted the same experiment on an icosehedral mesh of the unit sphere with = / = 1 



(following the "/-sphere" experiment of Thuburn et al. (2009)) and obtained the same result. 



4-2. Kelvin waves in a circular basin 

Coastal Kelvin waves provide a challenging test since they propagate at the gravity wave speed 
along the coast but are geostrophically balanced in the direction normal to the coast. We used 
the Kelvin wave initial condition for a circular basin with unit dimensionless radius as proposed in 



Ham et al. (2007), with Ro = 0.1 and Fr = 1. We integrated the equations until 10 dimensionless 
time units with a time step size At = 0.01. 

The mesh used for the Kelvin wave calculation is shown in Figure [2j Some snapshots of the 
numerical solution are shown in Figure [3j There are no spurious gravity waves observed, which 
means that the BDFMl discretisation is maintaining geostrophic balance in the normal direction 
as well as the Kelvin wave structure. 
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Figure 2: Mesh used for the Kelvin wave tests. 




Figure 3: Snapshots of the free surface elevation for the circular Kelvin wave testcase obtained at times t = 
0, 2500000, 5000000. The numerical scheme maintains the geostrophic balance in the normal direction, as indicated 
by the lack of radiated inertia-gravity waves. 
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Figure 4: Plot of errors from the Rossby convergence test with Rossby number Ro = le — 3 and timestep size 
At = 0.007996. The comparison is made after time tt/{1 + 87r^)/2 after which time the wave has travelled halfway 
around the domain. For large Ax we observe third-order convergence in both I2 and loo norms; for smaller Ax the 
error is dominated by either the timestepping error or the O(Ro^) truncation error in the small Rossby number 
expansion. 



4-3. Rossby waves 

To verify the convergence of the method we compared against the Rossby wave solution with 
streamfunction 

ip{x, y, t) = sin(27ra;) sin {2n {y + 7^)) , 7 



27r 



1 + 87r2' 

in a square domain with nondimensional length 1, with nondimensional wave propagation speed 
c = Ro^, and non-dimensional Coriolis parameter 



1 + Rot/ 
Ro ' 



and periodic boundary conditions in the x-direction. This is an exact solution of the Rossby 
wave equation, but is only an asymptotic limit solution of the linearised rotating shallow-water 
equations as Ro — )■ 0, with (9(Ro^) error. This means that for sufficiently small grid width and 
time step size we expect the (9(Ro^) error to dominate. The numerical solution was initialised 



from this streamfunction following the balanced initialisation approach described in Section AA_ 
A plot of the error is shown in figure El We observe 0{/A.x^ 



convergence until the error saturates 
because of the finite Rossby number. We attribute this third order convergence to the fact that 



in Section |2.8| the discrete Rossby wave equation was shown to be equal to usual continuous 
finite element discretisation of the Rossby wave equation using the space E, plus a perturbation. 
Since E contains all of the continuous piecewise quadratic functions, we would expect third-order 
convergence provided that the perturbation converges to zero sufficiently fast (although we do not 
currently have any estimates for the convergence of the perturbation). 
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Figure 5: Rossby waves on the "/3-tube" initialised from a streamfunction on a cylinder with a coarse un- 
structured triangle mesh. Colour plots of the free surface elevation are plotted at non-dimensional times 
t = 0.79957,19.9892,39.9784,59.9686,79.9568 from left to right. No unbalanced motions are visible from the 
plot. 



To demonstrate the performance of the numerical scheme on arbitrary manifolds we constructed 
an unstructured mesh of a cylinder with unit dimensionless radius and dimensionless height equal 
to 2. The Coriolis parameter was set to / = (1 -|- Ro2;)/Ro and other parameters were kept the 
same as the planar Rossby wave tests. We call this configuration the "/3"-tube since it corresponds 
to a /3-plane that has been wrapped into a cylinder. Some plots of the numerical integration of 
this test case are provided in Figure [5j no unbalanced motions are visible from the plots. 

4-4- Solid rotation on the sphere 

To investigate the grid imprinting caused by the finite element scheme, we integrated the 
linear rotating shallow-water equations on the sphere with initial condition obtained from the 
streamfunction ip = —uq cos9, where 9 is the latitude, uq = 27ri?/(12 days), and R = 6.37122 x 10^ 
is the radius of the sphere. The rotation rate \Q\ was 1/(1 day), and g = 9.8. This solution is a 
steady state solution of the linear equations with varying / because of the cylindrical symmetry; 
in general we do not expect numerical discretisations which break this symmetry to preserve the 
steady state. 

In our experiment, we used a level 4 icosahedral mesh (each icosahedron edge being subdivided 
into 8) of the sphere. The velocity and free surface elevation were initialised according to the 



procedure described in Section [4T| To measure the deviation from a steady state, the free surface 
elevation after 10 days of simulation with a timestep of 3600s was subtracted from the initial 
condition. Remarkably, the errors were almost indistinguishable from round-off error. It turns 
out that this is because of the mapping used between functions on the sphere, and functions on 
the icosahedral mesh with flat triangular elements used for the numerical integration. The flnite 
element streamfunction was initialised according to ip^ = ip o (p^ where is the mapping given 
as follows: 

^2 _ ^2x1/2 /^2_ ^2x1/2 
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Figure 6: Plots showing the exact steady numerical solution obtained using the balanced initialisation procedure. 
Top Left: The free surface elevation field. Top Right: The velocity field, plotted by evaluating the finite element 
field at vertices and edge midpoints of each triangle. Since only the normal components are continuous, there 
are multiple values of these vectors corresponding to the different elements that share those vertices/midpoints. 
Bottom: Close-up of the velocity vectors near the equator. 



This mapping preserves the value of z, and rescales x and y onto the sphere. Hence, we obtain 
ip^ = z, which can be represented exactly in the streamfunction space E. The same mapping is 
also applied to the finite element representation of the Coriolis parameter /, and we obtain 
f'^ = 2\Q\z which can also be represented exactly. Following the balanced initialisation procedure, 
the finite element free surface elevation field i]^ is obtained by projecting the mapping rj o into 
the pressure space V, where t] is the continuous balanced free surface elevation. Substitution into 
the velocity equation gives 

[ w^-u''dV 
dtJn 

[definition of f^, ip'^ and r]^] 
[integration by parts] 



where the second step follows since V ■ G V and so we can use the fact that t]^ is a finite 
element projection of rj in V, and where in the last step integration by parts was possible since ry 
is continuous and has continuous normal components. 



Jn Jn 

I fw^ - Vi^dV + [ V-io\d1/, 
Jn Jn 

[ - V I - c^ri] dV = 0, 
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5. Summary and outlook 



In this paper we described some properties of applying finite element spaces satisfying the div 
and curl embedding properties, applied to the rotating linear shallow water equations, in order 
to illustrate their possible suitability for numerical weather prediction on quasi-uniform grids. In 
this context, these methods can be thought of as more flexible extensions of the mimetic C-grid 
finite difference method that is currently used in many dynamical cores. This extra flexibility 
means that non-orthogonal grids and grids with rapid changes of mesh resolution can be used, 
and the ratio of pressure and velocity degrees of freedom can be adjusted to avoid spurious mode 
branches. We showed that spurious inertia-gravity mode branches will exist if dim(ii^) < dim(l^) 
and spurious Rossby mode branches will exist if dim(V) > dim(£'). The discrete Helmholtz 
decomposition implies that dim{E) = dim(S') — dim(V^) + 2 — dim(if) where H is the space of 
harmonic velocity fields on the chosen domain. This motivates the search for finite element spaces 
with dim(S') = 2 dim(V^) that can be used on pseudo-uniform grids on the sphere. In Section [s] 
we gave two low-order examples: the modified RTO-QO element pair for the non-orthogonal cubed 
sphere, and the BDFMI-PIdg element pair for triangles, the latter of which was illustrated with 
some numerical examples in Section |4} 

In future work, we shall aim to benchmark the augmented mixed element pair in the context 
of numerical weather prediction and ocean modelling. One of the benefits of this pair is that 
discontinuous Galerkin methods can be used for the nonlinear continuity equation for the density. 
These methods are locally conservative, have minimal dispersion and diffusion errors, and can 



be made TVB in combination with appropriate slope limiters as described in Cockburn and Shu 



(2001). Furthermore, as described in (White et al. , 2008), if one wishes to have tracer transport 



that is both conservative and consistent, it is necessary use the pressure space for tracers too. This 
means that tracer transport can (must) also use the discontinuous Galerkin method. 

Finally, we note that although the BDFM(k)-Pk£)G' finite element spaces do not have a 2:1 
ratio of velocity DOFs to pressure DOFs, there does exists a family of higher-order versions of 
the BDFMl element pair with a 2:1 ratio, obtained by appropriately augmenting the BDM(fc) 
spaces (with k > odd) with higher-order components that vanish on element boundaries. This 
does not work out so neatly for A; > 1 since it is also necessary to augment the P{k) space for 
pressure, to obtain stable element pairs with twice as many velocity DOF as pressure DOF per 
triangle. In future work, we shall investigate these higher-order element pairs, as well as extensions 
to tetrahedra in three-dimensions that can be used in unstructured mesh ocean modelling. 
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